% Process the eej strength derived from satellite data
% and compare it with ground deribed eej indices
%start date 30 OCT 2005
%latest date 21 JUNE 2006
%The latest changes were done for Stefan's new set of data sent on 20-6-6
load c:\manoj\projects\eej\EEJ_Stefan data;

% Note added on 26 June, 2006
% With the data from AAE-QSB (minute, 2001-2004), a total of 275 pairs were
% available in all the bins the correlations (without SQ & LT corr and with
% robust fitting) were 0.0785    0.6533    0.6061    0.7038    0.7021    0.6409    0.6088
% 0.7769   -0.3808
% The increase in the central bin with good significance (0.702) is evident


%Control Center
bin10 = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
obsdata = 7;  %ETT-HYB = 1,ETT-ABG = 2;TIR-HYB = 3;TIR-ABG = 4; HUA-FUQ 5; MBO-GUI = 6; 
%AAE-QSB 7; GUA-CBI 8; YAP-BIK 9; PND-HYB 10; ETT-PND 11; TIR - PND 12;
%AAE-ELT 13
obslt = [9.9 13.1]; %NOTE the LT are represented as for eg 8.99999 or 9.001 the filter
incrlat = 10;
mindate = tmjd(2001,01,01,0,0,0,1); % starting date
maxdate = tmjd(2004,12,31,0,0,0,1)+1; % ending date

%load c:\manoj\projects\eej\guacm4.mat index sq_est

switch obsdata,
    case 1,%ETT-HYB
        obslon = 82.5; %Oservatory Longitude of ETT is given as 82.5 as the IST is referenced to this longitude
        st = 12.8; % The bins are still centered around ETT 
        en = 132.8;

    case 2,%ETT-ABG
        obslon = 82.5;
        st = 12.8;
        en = 132.8;

    case 3,%TIR-HYB
        obslon = 82.5;
        st = 12.8;
        en = 132.8;

    case 4,%TIR-ABG
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
        plotoptions = 'k.-';

    case 5,%HUA-FUQ
        obslon = 285;
%        st = -142;
%        en = -8;
        st  = -139.5
        en = -10.5
        
    case 6,%MBO-GUI
        obslon = 343;
        st = -81.5;
        en = 47.5;
    case 7, %AAE-QSB
        obslon = 37.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
    case 8,%GUA-CBI
        obslon = 143;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
    case 9,%YAP=BIK
        obs_index = 2;
        obslon = 138.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
   case 10,%PND-HYB
        obslon = 82.5; %79 degrees 55 minutes
%         st = obslon - (incrlat*6+incrlat/2);
%         en = obslon + (incrlat*6+incrlat/2);
        st = 15; % The bins are still centered around ETT 
        en = 135;
   case 11,%ETT-PND
        obslon = 82.5; %Oservatory Longitude of ETT is given as 82.5 as the IST is referenced to this longitude
        st = 12.8; % The bins are still centered around ETT 
        en = 132.8;
   case 12,%TIR-PND
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
   case 13, %AAE-ELT
        obslon = 37.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
     otherwise,
        display('Error');
        %break;
        return;

end;


%INFORMATION
%     The EEJ parameters in the file EEJ.TXT (version JUNE 2006)
% 1.  FDAY
% 2.  geographic longitude of the magnetic dip equator
% 3.  geographic latitude of the magnetic dip equator
% 4.  local_time at the magnetic dip equator
% 5.  max current (old definition, assuming that we can estimate the zero line)
% 6.  av. current difference to north and south troughs 
% 7.  current difference to south trough 
% 8.  current difference to north trough 
% 9.  magnetic field strength at 108 km altitude
% 10. Kp
% 11. F10.7

%----------------------------


L = data(:,2) > st & data(:,2) < en+incrlat & data(:,1) < maxdate  &  data(:,1) > mindate...
    & data(:,12) >= obslt(1) & data(:,12) <= obslt(2) & data(:,10) <= 2.0 & data(:,5) >= 0.03;

load  c:\manoj\projects\eej\AAEQSBMIN20012004;
% a1 = repmat(aaehmean,[60*24,1]);
% a1 = reshape(a1,[1,60*24*1461]);
% a1(end-1439:end) = [];
% mjdax(1:1440) = [];
% aaeh(1:1440) = [];
% aaedeltah = aaeh-a1;
% 
% a1 = repmat(qsbhmean,[60*24,1]);
% a1 = reshape(a1,[1,60*24*1461]);
% a1(end-1439:end) = [];
% qsbh(1:1440) = [];
% qsbdeltah = qsbh-a1;




%The file contains variables (minute means of h)aaedeltah, qsbdeltah etc rom 2001-2004

sat_data = data(L,:);
sat_data(1:6,:) = [];%This is because as aae and qsb needs previous night's mean to remove
%the bias day 1.1.2001 need to be reomved from satellite and obs data set

for i = 1:length(sat_data),
    L = mjdax>=sat_data(i,1);
    kk(i) = sum(L);
end;

d=length(mjdax);

%Finding running mean over an hour for each pass
for i = 1:length(sat_data),
    eej_data_mean(i) = nanmean(aaedeltah(d-kk(i)-30:d-kk(i)+30)-...
        qsbdeltah(d-kk(i)-30:d-kk(i)+30));%running mean cetered on satellite pass
    eej_data_mom = aaedeltah(d-kk(i))-qsbdeltah(d-kk(i));%momentary ground data, w.r.t satellite pass
end;

% Process the eej strength derived from satellite data
% and compare it with ground deribed eej indices
%start date 30 OCT 2005
%latest date 29 MARCH 2006



%Control Center

x1data = bin10;
sqcor = 0; %SQ correction
ltcor = 1; %LT correction
robft = 1; %Robust outlier elimination
plotopt = 'k-';
plotcol = 'k';
plotchar = 'c';
ncoeff = 20;%same for hourly means 20 for minute means

%load c:\manoj\projects\eej\guacm4.mat index sq_est


        data1(:,1) = mjdax(d-kk);
        data1(:,2) = eej_data_mean;
        data1(isnan(data1(:,2)),2) = 9999; 
        data = sat_data;
        obslon = 37.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'k*-';

%----------------------------

ndd = zeros([1,25]);
nloop = 1;
lto = [];
ltc = [];
xdata = [];
nbindata = [];
ceff = [];
cerr = [];
nbin = 1;
obs_index = 2;
eej_index = 6;
nbin=1;
x1data = bin10;
clear ceff cerr;
nloop=1;
ndd = [];


for ki = st:incrlat:en,
L = data(:,2) > ki & data(:,2) < ki+incrlat ; %The last filter is to select only data before Dec 31, 2002
IndiaE = data(L,:);
ndata = 1;
ind1 = [];
ind2 = [];
ltt1 = [];
ltt2 = [];
nnodata=1;
for i = 1:length(IndiaE(:,1)),
    ind = findnearest(IndiaE(i,1),data1(:,1));
    [lt1,dummy1] = champ_lt(IndiaE(i,1),IndiaE(i,2),obslon);
    [dummy2,lt2] = champ_lt(data1(ind,1),IndiaE(i,2),obslon); % This is because observatory FDAY ~= CHAMP FDAY
    ltt1(i) = lt1;
    ltt2(i) = lt2;
    ltt3(i) = dummy1;
    ltt4(i) = dummy2;
    obst(i) = data1(ind,1);
    chpt(i) = IndiaE(i,1);
    
         if lt2 >= obslt(1) & lt2 <= obslt(2) & abs(data1(ind,obs_index)) < 300 & data1(ind,1) < max(data(:,1)),% & kp(ind) <= 2.5, %The last fliter take care the nodata (9999) points
            ind1(ndata) = ind;
%             f_ch(ndata) = polyval(p,lt1); %Getting the expexted EEJ strength from poly coefficients for champ LT
%             f_ob(ndata) = polyval(p,lt2); %Getting the expexted EEJ strength from poly coefficients for obs LT
            f_ch(ndata) = interp1(1:24,k,lt1,'spline'); %Found spline interpolation working better 18.1.06 !
            f_ob(ndata) = interp1(1:24,k,lt2,'spline');
%             fprintf('Lon = %5.1f i = %d f_ch = %f f_ob = %f lt1 = %f lt2 = %f kp = %3.1f\n', ki,i,f_ch(ndata),f_ob(ndata),lt1,lt2,data1(ind,end));
            ndata = ndata+1;
        else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
        end;
end;

if length(ind2) > 0,
    IndiaE(ind2,:) = [];
    obst(ind2) = [];
    chpt(ind2) = [];
end;

if length(ind1) > 1,
    obsE = data1(ind1,:);
    X = IndiaE(:,eej_index);
    if ltcor == 1,
        Y = obsE(:,obs_index).*(f_ch./f_ob)';
    elseif ltcor == 0,
        Y = obsE(:,obs_index)';
    end;
    
    if robft == 1,    
%    Robust outlier rejection  -- START
     if length(X) > 3,
     [rs,rstat] = robustfit(X,Y);
     residuals = abs(rstat.resid);
     LLL = residuals == max(residuals); %Just remove the most noisy point
     X(LLL) = [];
     Y(LLL) = [];
     end;
    % Robust outlier rejection  --- END
    end;

    [cof,err] = corrcoef(X,Y);  % correlation & Correction for observatory data 
    ceff(nloop) = cof(2,1);
    cerr(nloop) = err(2,1);
    fprintf('%4.3f ',cof(2,1));
    xdata(nloop) = ki+incrlat/2;
    ndd(nloop) = ndata-1;
    lto(nloop) = min(ltt2);
    ltc(nloop) = max(ltt2);
    nloop = nloop +1;

    %----------- PLOT SCATTER
else,
    x1data(nbin) = NaN;
     
end;
nbin = nbin + 1;
clear ind f_ch f_ob ind1 ind2 ndata nnodata IndiaE obsE ltt1 ltt2 ltt3 ltt4 obst chpt X Y;
end;

x1data(isnan(x1data)) = [];
fprintf('\n');
figure(1);
hold on;
errorbar(x1data,ceff,cerr/2,plotopt,'LineWidth',2);
axis([-80 80 -.6 1.4]);
grid on;
hold on;
